Guide to the effect of felzartamab anti-CD38 treatment on the molecular phenotype of antibody-mediated rejection in kidney transplant biopsies.

0.1 Foreword

This guide is meant to assist the interpretation of analyses carried out within “Diebold and Gauthier et al. Effect of felzartamab anti-CD38 treatment on the molecular phenotype of antibody-mediated rejection in kidney transplant biopsies. under review. Nature Medicine”. This manusript is currently under peer-review and has not yet been accepted for publication. This guide is a companion to the GitHub repository (https://github.com/TSI-PTG/CD38-effect-of-treatment)

1 Analysing the effect of treatment on molecular scores using aligned-rank transformed analysis of variance

1.1 INTRODUCTION

This analysis assesses how felzartamab therapy affects molecular scores measured with gene expression data from biopsies from kidney allografts taken from patients with antibody-mediated rejection. The assumptions of normality and heterogeneity of variance differ across molecular scores, thus we avail to aligned-rank transformation to carry out a non-parametric ANOVA. This is powerful in this context because it allows us to test for the interactive effect of treatment (i.e., placebo vs felzartamab) and time (i.e., change in scores from biopsy to biopsy). Aligned-rank transformation also allows for us to specify a mixed-model to handle repeated measurements (i.e., biopsies) across patients.

References:

• Wobbrock J, F.L., Gergle D, Higgins J The Aligned Rank Transform for Nonparametric Factorial Analyses Using Only ANOVA Procedures. In Proceedings of the ACM Conference on Human Factors in Computing Systems (CHI ’11) 143–146. doi:10.1145/1978942.1978963, https://depts.washington.edu/acelab/proj/art/.(2011).

• Kay M, E.L., Higgins J, Wobbrock J. ARTool: Aligned Rank Transform for Nonparametric Factorial ANOVAs. doi:10.5281/zenodo.594511, R package version 0.11.1, https://github.com/mjskay/ARTool. (2021).

• Elkin, L.A., Kay, M., Higgins, J.J. & Wobbrock, J.O. An Aligned Rank Transform Procedure for Multifactor Contrast Tests. in The 34th Annual ACM Symposium on User Interface Software and Technology 754–768 (Association for Computing Machinery, Virtual Event, USA, 2021).

1.2 DATA WRANGLING

First we need to load the necessary R packages and data.

# HOUSEKEEPING ####
# Load necessary libraries
library(tidyverse)

## -- Attaching core tidyverse packages -------------------- tidyverse 2.0.0 --
## v dplyr     1.1.4     v readr     2.1.5
## v forcats   1.0.0     v stringr   1.5.1
## v ggplot2   3.5.1     v tibble    3.2.1
## v lubridate 1.9.3     v tidyr     1.3.1
## v purrr     1.0.2     
## -- Conflicts -------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(broom)
library(rstatix)

## 
## Attaching package: 'rstatix'
## 
## The following object is masked from 'package:stats':
## 
##     filter

library(flextable)

## 
## Attaching package: 'flextable'
## 
## The following object is masked from 'package:purrr':
## 
##     compose

library(officer)
library(emmeans)

## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'

library(multcomp)

## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:rstatix':
## 
##     select
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 
## Attaching package: 'TH.data'
## 
## The following object is masked from 'package:MASS':
## 
##     geyser

library(ARTool)
library(rcompanion)
library(Biobase)

## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## 
## The following object is masked from 'package:flextable':
## 
##     width
## 
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## 
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## 
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter,
##     Find, get, grep, grepl, intersect, is.unsorted, lapply, Map,
##     mapply, match, mget, order, paste, pmax, pmax.int, pmin,
##     pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
##     setdiff, table, tapply, union, unique, unsplit, which.max,
##     which.min
## 
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

library(ggpubr) 

## 
## Attaching package: 'ggpubr'
## 
## The following objects are masked from 'package:flextable':
## 
##     border, font, rotate

library(patchwork) 

## 
## Attaching package: 'patchwork'
## 
## The following object is masked from 'package:MASS':
## 
##     area

library(ggprism) 

# Custom operators
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0

# Suppress dplyr reframe info
options(dplyr.reframe.inform = FALSE)

# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/results/flextable_artANOVA.RData")
# load plot data
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_artANOVA.RData")

# Seed for reproducibility
set.seed(42)

Now we define the categories of molecular scores for easier organization downstream.

# Feature categories
vars_cfDNA <- c("cfDNA_cpml")
vars_abmr <- c("ABMRpm", "ggt0", "ptcgt0", "DSAST", "NKB")
vars_tcmr <- c("TCMRt", "tgt1", "igt1", "QCAT", "TCB")
vars_injury <- c("IRRAT30", "IRITD3", "IRITD5", "cigt1", "ctgt1")
vars_parenchyma <- c("KT1", "KT2")
vars_macrophage <- c("AMAT1", "QCMAT")

# DEFINE VARIABLES TO ASSESS ####
vars <- c(vars_cfDNA, vars_abmr, vars_tcmr, vars_macrophage, vars_injury, vars_parenchyma)

Now we process the data so we can iteratively carry out the analyses on each score.

# DEFINE THE DATA ####
data <- set %>%
    pData()  %>% 
    dplyr::select(Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, all_of(vars)) %>%
    dplyr::filter(Patient %nin% c(15, 18)) %>%
    left_join(., summarise(., sample_pairs = n(), .by = Felzartamab_Followup), by = "Felzartamab_Followup") %>%
    relocate(sample_pairs, .after = Felzartamab_Followup) %>%
    arrange(Felzartamab, Patient, Group)


# WRANGLE THE PHENOTYPE DATA ####
data_00 <- data %>%
    expand_grid(category = c("cfDNA", "ABMR", "TCMR", "macrophage", "injury", "parenchyma")) %>%
    nest(.by = category) %>%
    mutate(
        features = map(
            category,
            function(category) {
                if (category == "cfDNA") {
                    vars_cfDNA
                } else if (category == "ABMR") {
                    vars_abmr
                } else if (category == "TCMR") {
                    vars_tcmr
                } else if (category == "macrophage") {
                    vars_macrophage
                } else if (category == "injury") {
                    vars_injury
                } else if (category == "parenchyma") {
                    vars_parenchyma
                } 
            }
        ),
        data = pmap(
            list(features, data),
            function(features, data) {
                data %>%
                    dplyr::select(
                        Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, sample_pairs,
                        all_of(features)
                    )
            }
        )
    )


# WRANGLE THE DATA FOR UNIVARIATE TESTS ####
data_01 <- data_00 %>%
    mutate(
        data_univariate = pmap(
            list(data, features),
            function(data, features) {
                data %>%
                    pivot_longer(
                        cols = all_of(features),
                        names_to = "variable",
                        values_to = "value"
                    ) %>%
                    nest(.by = variable) %>%
                    mutate(
                        variable = variable %>%
                            factor(
                                levels = vars
                            ),
                        annotation = case_when(
                            variable %in% vars_cfDNA ~ "cfDNA",
                            variable %in% vars_tcmr ~ "TCMR-related",
                            variable %in% vars_abmr ~ "ABMR-related",
                            variable %in% vars_macrophage ~ "macrophage-related",
                            variable %in% vars_injury ~ "injury-related",
                            variable %in% vars_parenchyma ~ "parenchyma-related",
                            TRUE ~ " "
                        ) %>%
                            factor(
                                levels = c(
                                    "cfDNA",
                                    "ABMR-related",
                                    "TCMR-related",
                                    "macrophage-related",
                                    "injury-related",
                                    "parenchyma-related",
                                    "archetypes",
                                    "rejection PC",
                                    "injury PC"
                                )
                            ),
                        score = case_when(
                            variable == "cfDNA_cpml" ~ "Donor-derived cell-free DNA (dd-cfDNA, cp/mL)",
                            variable == "TCMRt" ~ "TCMR classifier (TCMRProb)",
                            variable == "TCB" ~ "T cell burden (TCB)",
                            variable == "tgt1" ~ "Tubulitis classifier (t>1Prob)",
                            variable == "igt1" ~ "Interstitial infiltrate classifier (i>1Prob)",
                            variable == "QCAT" ~ "Cytotoxic T cell-associated (QCAT)",
                            variable == "ABMRpm" ~ "ABMR classifier (ABMRProb)",
                            variable == "DSAST" ~ "DSA-selective transcripts (DSAST)",
                            variable == "NKB" ~ "NK cell burden (NKB)",
                            variable == "ggt0" ~ "Glomerulitis classifier (g>0Prob)",
                            variable == "ptcgt0" ~ "Peritubular capillaritis classifier (ptc>0Prob)",
                            variable == "AMAT1" ~ "Alternatively activated macrophage (AMAT1)",
                            variable == "QCMAT" ~ "Constitutive macrophage (QCMAT)",
                            variable == "IRITD3" ~ "Injury-repair induced, day 3 (IRITD3)",
                            variable == "IRITD5" ~ "Injury-repair induced, day 5 (IRITD5)",
                            variable == "IRRAT30" ~ "Injury-repair associated (IRRAT30)",
                            variable == "cigt1" ~ "Fibrosis classifier (ci>1Prob)",
                            variable == "ctgt1" ~ "Atrophy classifier (ct>1Prob)",
                            variable == "KT1" ~ "Kidney parenchymal (KT1)",
                            variable == "KT2" ~ "Kidney parenchymal - no solute carriers (KT2)"
                        ), .before = 1
                    ) %>%
                    arrange(annotation, variable)
            }
        )
    ) %>%
    dplyr::select(category, data_univariate) %>%
    unnest(data_univariate) %>%
    mutate(n_cat = score %>% unique() %>% length(), .by = category, .after = variable)

This is what the data look like after we have organized them. Each row of the dataframe contains the data for each score, the category it belongs to, its category annotation, its long-format name, and the name of the variable in the data.

data_01 %>% print(n="all")

## # A tibble: 20 x 6
##    category   annotation         score              variable n_cat data    
##    <chr>      <fct>              <chr>              <fct>    <int> <list>  
##  1 cfDNA      cfDNA              Donor-derived cel~ cfDNA_c~     1 <tibble>
##  2 ABMR       ABMR-related       ABMR classifier (~ ABMRpm       5 <tibble>
##  3 ABMR       ABMR-related       Glomerulitis clas~ ggt0         5 <tibble>
##  4 ABMR       ABMR-related       Peritubular capil~ ptcgt0       5 <tibble>
##  5 ABMR       ABMR-related       DSA-selective tra~ DSAST        5 <tibble>
##  6 ABMR       ABMR-related       NK cell burden (N~ NKB          5 <tibble>
##  7 TCMR       TCMR-related       TCMR classifier (~ TCMRt        5 <tibble>
##  8 TCMR       TCMR-related       Tubulitis classif~ tgt1         5 <tibble>
##  9 TCMR       TCMR-related       Interstitial infi~ igt1         5 <tibble>
## 10 TCMR       TCMR-related       Cytotoxic T cell-~ QCAT         5 <tibble>
## 11 TCMR       TCMR-related       T cell burden (TC~ TCB          5 <tibble>
## 12 macrophage macrophage-related Alternatively act~ AMAT1        2 <tibble>
## 13 macrophage macrophage-related Constitutive macr~ QCMAT        2 <tibble>
## 14 injury     injury-related     Injury-repair ass~ IRRAT30      5 <tibble>
## 15 injury     injury-related     Injury-repair ind~ IRITD3       5 <tibble>
## 16 injury     injury-related     Injury-repair ind~ IRITD5       5 <tibble>
## 17 injury     injury-related     Fibrosis classifi~ cigt1        5 <tibble>
## 18 injury     injury-related     Atrophy classifie~ ctgt1        5 <tibble>
## 19 parenchyma parenchyma-related Kidney parenchyma~ KT1          2 <tibble>
## 20 parenchyma parenchyma-related Kidney parenchyma~ KT2          2 <tibble>

1.3 MEDIAN SCORES AND MEDIAN EFFECT OF TREATMENT

At this stage we can summarize the median values for each score by treatment and follow-up.

(Note that the median scores are used for general interpretation of treatment effects because aligned-ranks, the data being assessed for inference regarding the effects of treatment, are not easily interpretable.)

# UNIVARIATE MEDIANS ####
data_02 <- data_01 %>%
    mutate(
        medians = map(
            data,
            function(data) {
                data %>%
                    reframe(
                        median = value %>% median(),
                        IQR = value %>% IQR(),
                        sample_pairs = n(),
                        .by = c(Followup, Felzartamab, Felzartamab_Followup)
                    ) %>%
                    arrange(Felzartamab_Followup)
            }
        ),
        medians_delta = map(
            medians,
            function(medians) {
                medians %>%
                    reframe(
                        median_delta = combn(median, 2, diff) %>% as.numeric(),
                        IQR_delta = combn(IQR, 2, mean) %>% as.numeric(),
                        sample_pairs = sample_pairs,
                        .by = c(Felzartamab)
                    ) %>%
                    mutate(
                        Followup_pairwise = rep(c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52"), 2) %>%
                            factor(levels = c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52")),
                        .before = sample_pairs
                    ) %>%
                    mutate(
                        median_delta_delta = combn(median_delta, 2, diff) %>% as.numeric(),
                        IQR_delta_delta = combn(IQR_delta, 2, mean) %>% as.numeric(),
                        .by = c(Followup_pairwise),
                        .before = sample_pairs
                    )
            }
        )
    )

Let’s take a look at what the median summaries look like for the ABMRprob score.

data_02 %>%
    dplyr::filter(variable == "ABMRpm") %>%
    pull(medians) %>%
    pluck(1) %>%
    flextable::flextable()

Followup	Felzartamab	Felzartamab_Followup	median	IQR	sample_pairs
Baseline	Placebo	Baseline_Placebo	0.6695	0.35525	10
Week24	Placebo	Week24_Placebo	0.7740	0.39025	10
Week52	Placebo	Week52_Placebo	0.5185	0.49825	10
Baseline	Felzartamab	Baseline_Felzartamab	0.7360	0.30775	10
Week24	Felzartamab	Week24_Felzartamab	0.1675	0.35700	10
Week52	Felzartamab	Week52_Felzartamab	0.7025	0.53875	10

Now we can look at the change in median ABMRprob scores in the placebo and felzartamab arms, as well as the difference in those medians (i.e., ΔΔ median - this is our summarized treatment effect).

data_02 %>%
    dplyr::filter(variable == "ABMRpm") %>%
    pull(medians_delta) %>%
    pluck(1) %>%
    flextable::flextable()

Felzartamab	median_delta	IQR_delta	Followup_pairwise	median_delta_delta	IQR_delta_delta	sample_pairs
Placebo	0.1045	0.372750	Baseline - Week24	-0.6730	0.3525625	10
Placebo	-0.1510	0.426750	Baseline - Week52	0.1175	0.4250000	10
Placebo	-0.2555	0.444250	Week24 - Week52	0.7905	0.4460625	10
Felzartamab	-0.5685	0.332375	Baseline - Week24	-0.6730	0.3525625	10
Felzartamab	-0.0335	0.423250	Baseline - Week52	0.1175	0.4250000	10
Felzartamab	0.5350	0.447875	Week24 - Week52	0.7905	0.4460625	10

1.4 ALIGNED-RANK TRANSFORM ANOVA

Now we carry out the analysis. This includes testing the specific contrasts for the interaction of treatment and time (i.e., our effect of treatment).

# UNIVARIATE NONPARAMETRIC TESTS ####
artANOVA <- data_02 %>%
    mutate(
        art = map(data, art, formula = value ~ Followup * Felzartamab + (1 | Patient)),
        art_aov = map(art, anova, type = "II"),
        art_aov_tidy = map(art_aov, tidy) %>% suppressWarnings(),
        art_aov_contrast = map(
            art,
            art.con,
            formula = "Followup:Felzartamab", adjust = "fdr", method = "pairwise", interaction = TRUE, response = "art"
        ),
        art_aov_contrast_tidy = map(art_aov_contrast, tidy),
        art_aov_contrast_cld = map(
            art_aov_contrast_tidy,
            function(art_aov_contrast_tidy) {
                art_aov_contrast_tidy %>%
                    as.data.frame() %>%
                    cldList(adj.p.value ~ Followup:Felzartamab, data = .)
            }
        )
    )

1.5 SUMMARIZE RESULTS FOR EACH SCORE

# CREATE FLEXTABLE SUMMARY OF ART MODELS ####
title_art <- paste("Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients")
artANOVA %>%
    dplyr::select(annotation, n_cat, score, art_aov) %>%
    unnest(everything()) %>%
    dplyr::rename(p.value = `Pr(>F)`) %>%
    mutate(FDR = p.value %>% p.adjust(method = "fdr"), .by = annotation) %>%
    dplyr::select(annotation, score, Term, `F`, p.value, FDR) %>%
    flextable::flextable() %>%
    flextable::add_header_row(values = rep(title_art, ncol_keys(.))) %>%
    flextable::merge_h(part = "header") %>%
    flextable::merge_v(j = 1:2) %>%
    flextable::fontsize(size = 8, part = "all") %>%
    flextable::align(align = "center", part = "all") %>%
    flextable::bg(bg = "white", part = "all") %>%
    flextable::bg(i = ~ FDR < 0.05 & Term == "Followup:Felzartamab", j = 3:6, bg = "#fbff00") %>%
    flextable::colformat_double(j = 2:3, digits = 2) %>%
    flextable::colformat_double(j = 4:ncol_keys(.), digits = 3) %>%
    flextable::border_remove() %>%
    flextable::bold(part = "header") %>%
    flextable::padding(padding = 0, part = "all") %>%
    flextable::border(border = officer::fp_border(), part = "all")  %>%
    flextable::autofit()

Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients
annotation	score	Term	F	p.value	FDR
cfDNA	Donor-derived cell-free DNA (dd-cfDNA, cp/mL)	Followup	4.418	0.019	0.029
		Felzartamab	1.598	0.222	0.222
		Followup:Felzartamab	10.396	0.000	0.001
ABMR-related	ABMR classifier (ABMRProb)	Followup	6.659	0.003	0.007
		Felzartamab	1.620	0.219	0.235
		Followup:Felzartamab	9.605	0.000	0.001
	Glomerulitis classifier (g>0Prob)	Followup	11.768	0.000	0.001
		Felzartamab	1.511	0.235	0.235
		Followup:Felzartamab	10.743	0.000	0.001
	Peritubular capillaritis classifier (ptc>0Prob)	Followup	14.188	0.000	0.000
		Felzartamab	3.178	0.092	0.125
		Followup:Felzartamab	7.749	0.002	0.004
	DSA-selective transcripts (DSAST)	Followup	10.508	0.000	0.001
		Felzartamab	1.747	0.203	0.234
		Followup:Felzartamab	4.179	0.023	0.035
	NK cell burden (NKB)	Followup	6.512	0.004	0.007
		Felzartamab	1.956	0.179	0.224
		Followup:Felzartamab	5.114	0.011	0.018
TCMR-related	TCMR classifier (TCMRProb)	Followup	3.555	0.039	0.165
		Felzartamab	1.417	0.249	0.468
		Followup:Felzartamab	1.124	0.336	0.504
	Tubulitis classifier (t>1Prob)	Followup	1.016	0.372	0.506
		Felzartamab	1.853	0.190	0.408
		Followup:Felzartamab	0.065	0.937	0.968
	Interstitial infiltrate classifier (i>1Prob)	Followup	3.414	0.044	0.165
		Felzartamab	0.997	0.331	0.504
		Followup:Felzartamab	0.032	0.968	0.968
	Cytotoxic T cell-associated (QCAT)	Followup	5.693	0.007	0.107
		Felzartamab	4.118	0.057	0.172
		Followup:Felzartamab	0.927	0.405	0.506
	T cell burden (TCB)	Followup	4.360	0.020	0.151
		Felzartamab	2.868	0.108	0.269
		Followup:Felzartamab	0.301	0.742	0.856
macrophage-related	Alternatively activated macrophage (AMAT1)	Followup	1.549	0.226	0.304
		Felzartamab	1.047	0.320	0.320
		Followup:Felzartamab	2.406	0.105	0.209
	Constitutive macrophage (QCMAT)	Followup	3.378	0.045	0.209
		Felzartamab	3.100	0.095	0.209
		Followup:Felzartamab	1.426	0.253	0.304
injury-related	Injury-repair associated (IRRAT30)	Followup	0.055	0.946	0.946
		Felzartamab	3.842	0.066	0.223
		Followup:Felzartamab	1.915	0.162	0.347
	Injury-repair induced, day 3 (IRITD3)	Followup	0.174	0.841	0.946
		Felzartamab	3.268	0.087	0.223
		Followup:Felzartamab	3.125	0.056	0.223
	Injury-repair induced, day 5 (IRITD5)	Followup	0.073	0.930	0.946
		Felzartamab	1.088	0.311	0.583
		Followup:Felzartamab	2.586	0.089	0.223
	Fibrosis classifier (ci>1Prob)	Followup	0.360	0.700	0.876
		Felzartamab	3.482	0.078	0.223
		Followup:Felzartamab	0.938	0.401	0.668
	Atrophy classifier (ct>1Prob)	Followup	0.568	0.572	0.780
		Felzartamab	3.297	0.086	0.223
		Followup:Felzartamab	0.602	0.553	0.780
parenchyma-related	Kidney parenchymal (KT1)	Followup	0.968	0.389	0.518
		Felzartamab	1.628	0.218	0.518
		Followup:Felzartamab	0.605	0.552	0.552
	Kidney parenchymal - no solute carriers (KT2)	Followup	0.859	0.432	0.518
		Felzartamab	0.911	0.352	0.518
		Followup:Felzartamab	1.038	0.365	0.518

1.6 SUMMARISE THE EFFECT OF TREATMENT FOR EACH SCORE ACROSS ALL TIME PERIODS

Here we present the effect of treatment (ΔΔ) for each score, along with the effect of time (Δ) for placebo and felzartmab arms.

flextable_pairwise

Table 1. Effect of felzartamab on molecular scores
Annotation	Score	Baseline - Week 24				Week 24 - Week 52				Baseline - Week 52
		Δ Placebo (N=10)	Δ Felzartamab (N=10)	ΔΔ	ΔΔ FDR	Δ Placebo (N=10)	Δ Felzartamab (N=10)	ΔΔ	ΔΔ FDR	Δ Placebo (N=10)	Δ Felzartamab (N=10)	ΔΔ	ΔΔ FDR
cfDNA	Donor-derived cell-free DNA (dd-cfDNA, cp/mL)	0.50 (76.12)	-55.00 (30)	-55.50 (53.06)	2e-04	-3.50 (69.12)	14.50 (25.5)	18.00 (47.31)	0.027	-3.00 (51)	-40.50 (45.75)	-37.50 (48.38)	0.045
ABMR-related	ABMR classifier (ABMRProb)	0.10 (0.37)	-0.57 (0.33)	-0.67 (0.35)	1e-03	-0.26 (0.44)	0.54 (0.45)	0.79 (0.45)	1e-03	-0.15 (0.43)	-0.03 (0.42)	0.12 (0.43)	0.709
	Glomerulitis classifier (g>0Prob)	-0.02 (0.23)	-0.31 (0.35)	-0.29 (0.29)	3e-03	-0.19 (0.34)	0.31 (0.41)	0.50 (0.37)	2e-04	-0.22 (0.3)	0.00 (0.38)	0.21 (0.34)	0.277
	Peritubular capillaritis classifier (ptc>0Prob)	-0.12 (0.25)	-0.55 (0.27)	-0.43 (0.26)	3e-03	-0.10 (0.37)	0.37 (0.31)	0.47 (0.34)	3e-03	-0.22 (0.3)	-0.18 (0.35)	0.04 (0.33)	0.977
	DSA-selective transcripts (DSAST)	-0.09 (0.22)	-0.37 (0.24)	-0.28 (0.23)	0.026	-0.16 (0.32)	0.16 (0.33)	0.32 (0.33)	0.026	-0.24 (0.33)	-0.20 (0.25)	0.04 (0.29)	0.988
	NK cell burden (NKB)	0.02 (0.22)	-0.78 (0.49)	-0.80 (0.36)	0.014	-0.19 (0.4)	0.68 (0.59)	0.87 (0.49)	0.014	-0.17 (0.46)	-0.10 (0.62)	0.07 (0.54)	0.972
TCMR-related	TCMR classifier (TCMRProb)	-0.01 (0.02)	0.00 (0)	0.01 (0.01)	0.358	0.00 (0.02)	0.00 (0)	0.00 (0.01)	0.856	-0.01 (0.02)	0.00 (0)	0.01 (0.01)	0.358
	Tubulitis classifier (t>1Prob)	-0.01 (0.06)	-0.02 (0.04)	-0.01 (0.05)	1.000	0.00 (0.04)	0.01 (0.02)	0.01 (0.03)	1.000	-0.01 (0.05)	-0.01 (0.03)	-0.01 (0.04)	1.000
	Interstitial infiltrate classifier (i>1Prob)	-0.01 (0.04)	-0.02 (0.03)	-0.01 (0.04)	0.964	0.01 (0.04)	0.01 (0.04)	0.00 (0.04)	0.964	0.00 (0.05)	-0.01 (0.06)	-0.01 (0.05)	0.964
	Cytotoxic T cell-associated (QCAT)	-0.20 (0.41)	-0.70 (0.46)	-0.50 (0.43)	0.425	0.10 (0.33)	0.54 (0.57)	0.44 (0.45)	0.425	-0.10 (0.46)	-0.16 (0.41)	-0.06 (0.44)	0.871
	T cell burden (TCB)	-0.44 (0.71)	-0.79 (0.56)	-0.36 (0.64)	0.825	0.28 (0.57)	0.34 (0.6)	0.06 (0.58)	0.905	-0.16 (0.88)	-0.45 (0.43)	-0.29 (0.65)	0.825
macrophage-related	Alternatively activated macrophage (AMAT1)	0.08 (0.2)	-0.27 (0.45)	-0.35 (0.33)	0.128	0.09 (0.23)	0.11 (0.47)	0.02 (0.35)	0.812	0.18 (0.19)	-0.16 (0.32)	-0.33 (0.25)	0.128
	Constitutive macrophage (QCMAT)	0.00 (0.15)	-0.24 (0.27)	-0.24 (0.21)	0.340	0.01 (0.12)	0.09 (0.27)	0.08 (0.2)	0.685	0.01 (0.11)	-0.15 (0.24)	-0.16 (0.18)	0.348
injury-related	Injury-repair associated (IRRAT30)	0.18 (0.52)	0.05 (0.4)	-0.13 (0.46)	0.353	0.05 (0.66)	-0.11 (0.72)	-0.16 (0.69)	0.353	0.23 (0.6)	-0.06 (0.48)	-0.29 (0.54)	0.175
	Injury-repair induced, day 3 (IRITD3)	0.03 (0.19)	0.04 (0.12)	0.01 (0.15)	0.945	0.03 (0.18)	-0.18 (0.18)	-0.20 (0.18)	0.060	0.06 (0.19)	-0.13 (0.14)	-0.19 (0.16)	0.060
	Injury-repair induced, day 5 (IRITD5)	-0.01 (0.17)	0.04 (0.14)	0.05 (0.15)	0.932	0.07 (0.2)	-0.17 (0.21)	-0.24 (0.21)	0.093	0.06 (0.18)	-0.13 (0.17)	-0.18 (0.18)	0.093
	Fibrosis classifier (ci>1Prob)	-0.02 (0.29)	0.13 (0.32)	0.15 (0.3)	0.434	0.06 (0.25)	-0.08 (0.53)	-0.14 (0.39)	0.844	0.04 (0.25)	0.05 (0.37)	0.01 (0.31)	0.434
	Atrophy classifier (ct>1Prob)	0.01 (0.33)	0.12 (0.31)	0.12 (0.32)	0.587	0.09 (0.27)	-0.06 (0.44)	-0.15 (0.36)	0.587	0.10 (0.32)	0.07 (0.41)	-0.03 (0.36)	0.587
parenchyma-related	Kidney parenchymal (KT1)	-0.05 (0.22)	-0.15 (0.23)	-0.09 (0.22)	0.731	-0.01 (0.11)	0.10 (0.22)	0.11 (0.17)	0.705	-0.06 (0.25)	-0.04 (0.18)	0.02 (0.22)	0.705
	Kidney parenchymal - no solute carriers (KT2)	-0.09 (0.39)	-0.16 (0.37)	-0.07 (0.38)	0.577	-0.09 (0.2)	0.14 (0.38)	0.23 (0.29)	0.484	-0.18 (0.42)	-0.02 (0.31)	0.16 (0.37)	0.577
Grey shading denotes ANOVA interactive effect FDR < 0.05 FDR correction was carried out within each annotation grouping

1.7 PLOT THE RESULTS

# EXTRACT LEGEND FOR PLOTS ####
panel_legend <- plots_artANOVA %>%
    dplyr::filter(variable == "AMAT1") %>%
    pull(plot_violin) %>%
    ggpubr::get_legend() %>%
    ggpubr::as_ggplot() +
    theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))


# MOLECULAR ABMR PANELS ####
panel_violin_abmr <- plots_artANOVA %>%
    dplyr::filter(category %in% c("ABMR")) %>%
    pull(plot_violin) %>%
    wrap_plots(nrow = 1, ncol = 5) &
    theme(
        axis.text = element_text(size = 10, colour = "black"),
        legend.position = "none",
        plot.background = element_rect(fill = "grey95", colour = " white")
    )

panel_violin_abmr <- panel_violin_abmr %>%
    ggarrange(
        labels = "B",
        font.label = list(size = 25, face = "bold"),
        legend = "none"
    ) %>%
    ggpubr::annotate_figure(
        top = text_grob("Effect of felzartamab on molecular ABMR activity scores",
            face = "bold.italic",
            size = 25,
            hjust = 1.27
        )
    )


# MOLECULAR TCMR PANELS ####
panel_violin_tcmr <- plots_artANOVA %>%
    dplyr::filter(category %in% c("TCMR")) %>%
    pull(plot_violin) %>%
    wrap_plots(nrow = 1, ncol = 5) &
    theme(
        axis.text = element_text(size = 10, colour = "black"),
        legend.position = "none",
    )

panel_violin_tcmr <- panel_violin_tcmr %>%
    ggarrange(
        labels = "C",
        font.label = list(size = 25, face = "bold"),
        legend = "none"
    ) %>%
    ggpubr::annotate_figure(
        top = text_grob("Effect of felzartamab on molecular TCMR activity scores",
            face = "bold.italic",
            size = 25,
            hjust = 1.27
        )
    )


# COMBINED VIOLIN PANELS ####
panels_violin <- ggarrange(
    panel_violin_abmr,
    panel_violin_tcmr,
    nrow = 2,
    heights = c(1, 1)
)

ggarrange(
    panel_legend,
    panels_violin,
    nrow = 2,
    heights = c(0.125, 1)
) 

2 Slope analysis of the effect of felzartamab on injury scores

2.1 INTRODUCTION

This analysis assesses the slopes for injury scores over time between placebo and felzartamab arms. We used a linear mixed-model to account for repeated measurements over time.

2.2 DATA WRANGLING

First we need to load the necessary R packages and data.

# HOUSEKEEPING ####
# CRAN packages
library(tidyverse)
library(lme4)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

library(lmerTest)

## 
## Attaching package: 'lmerTest'

## The following object is masked from 'package:lme4':
## 
##     lmer

## The following object is masked from 'package:stats':
## 
##     step

library(performance)
library(ggeffects)
library(flextable)
library(ggpubr) 
library(patchwork) 
library(ggprism) 
# Bioconductor libraries
library(Biobase) # BiocManager::install("Biobase")
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
source("C:/R/CD38-effect-of-treatment/natmed/code/functions/get_slope_function.r")
# load reference set
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
# load refression plots
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_regression.RData")
# load violin plots
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_artANOVA.RData")

# Seed for reproducibility
set.seed(42)


# DEFINE THE DATA ####
data <- set %>%
    pData() %>%
    dplyr::select(Center, Patient, Felzartamab, Group, IRRAT30, IRITD3, IRITD5) %>%
    dplyr::filter(Patient %nin% c(15, 18)) %>%
    mutate(
        time = case_when(
            Group == "Index" ~ 0,
            Group == "FU1" ~ 0.46,
            Group == "FU2" ~ 0.99
        ),
        linetype = "solid" %>% factor()
    )

Now we fit the linear mixed-models.

# FIT MODELS ####
fit_IRRAT30 <- lmer(IRRAT30 ~ Felzartamab * time + (time | Patient), data)
fit_IRITD3 <- lmer(IRITD3 ~ Felzartamab * time + (time | Patient), data)
fit_IRITD5 <- lmer(IRITD5 ~ Felzartamab * time + (time | Patient), data)

We check the validatity of the model assumptions for IRRAT30.

# test model assumptions
performance::check_model(fit_IRRAT30)

performance::check_normality(fit_IRRAT30)

## OK: residuals appear as normally distributed (p = 0.605).

We check the validatity of the model assumptions for IRITD3.

# test model assumptions
performance::check_model(fit_IRITD3)

performance::check_normality(fit_IRITD3)

## OK: residuals appear as normally distributed (p = 0.070).

We check the validatity of the model assumptions for IRITD5.

# test model assumptions
performance::check_model(fit_IRITD5)

performance::check_normality(fit_IRITD5)

## OK: residuals appear as normally distributed (p = 0.090).

2.3 ORGANIZE THE RESULTS

The next step is to extract and summarize the model results.

# EXTRACT THE SLOPES FROM THE MODEL FITS ####
# IRRAT30
# Create list object into which to place model results
slopes_IRRAT30 <- list()
# Acute slopes
slopes_IRRAT30[[1]] <- get_slope(
    .model_obj = fit_IRRAT30, .name = "Placebo",
    .contrasts = c(0, 0, 1, 0) #
)
slopes_IRRAT30[[2]] <- get_slope(
    .model_obj = fit_IRRAT30, .name = "Felzartamab",
    .contrasts = c(0, 0, 1, 1)
)
slopes_IRRAT30[[3]] <- get_slope(
    .model_obj = fit_IRRAT30, .name = "Intergroup Difference",
    .contrasts = c(0, 0, 0, 1)
)
model_IRRAT30 <- slopes_IRRAT30 %>%
    bind_rows() %>%
    mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
    dplyr::select(name, result, p_value, p_adjusted)

# IRITD3
# Create list object into which to place model results
slopes_IRITD3 <- list()
# Acute slopes
slopes_IRITD3[[1]] <- get_slope(
    .model_obj = fit_IRITD3, .name = "Placebo",
    .contrasts = c(0, 0, 1, 0) #
)
slopes_IRITD3[[2]] <- get_slope(
    .model_obj = fit_IRITD3, .name = "Felzartamab",
    .contrasts = c(0, 0, 1, 1)
)
slopes_IRITD3[[3]] <- get_slope(
    .model_obj = fit_IRITD3, .name = "Intergroup Difference",
    .contrasts = c(0, 0, 0, 1)
)
model_IRITD3 <- slopes_IRITD3 %>%
    bind_rows() %>%
    mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
    dplyr::select(name, result, p_value, p_adjusted)

# IRITD5
# Create list object into which to place model results
slopes_IRITD5 <- list()
# Acute slopes
slopes_IRITD5[[1]] <- get_slope(
    .model_obj = fit_IRITD5, .name = "Placebo",
    .contrasts = c(0, 0, 1, 0) #
)
slopes_IRITD5[[2]] <- get_slope(
    .model_obj = fit_IRITD5, .name = "Felzartamab",
    .contrasts = c(0, 0, 1, 1)
)
slopes_IRITD5[[3]] <- get_slope(
    .model_obj = fit_IRITD5, .name = "Intergroup Difference",
    .contrasts = c(0, 0, 0, 1)
)
model_IRITD5 <- slopes_IRITD5 %>%
    bind_rows() %>%
    mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
    dplyr::select(name, result, p_value, p_adjusted)

2.4 SUMMARIZE RESULTS

IRRAT30

size_font <- 10
model_IRRAT30 %>%
    dplyr::select(-p_value) %>%
    flextable::qflextable() %>%
    flextable::set_table_properties(layout = "autofit") %>%
    flextable::set_header_labels(
        name = "Group", result = "IRRAT30 Slope (95%CI)",
        "p_adjusted" = "FDR"
    ) %>%
    flextable::bold(i = NULL, part = "header") %>%
    flextable::fontsize(size = size_font, part = "all")

Group	IRRAT30 Slope (95%CI)	FDR
Placebo	0.32 (-0.04 to 0.67)	0.1133
Felzartamab	-0.29 (-0.64 to 0.07)	0.1133
Intergroup Difference	-0.60 (-1.10 to -0.10)	0.0546

IRITD3

model_IRITD3 %>%
    dplyr::select(-p_value) %>%
    flextable::qflextable() %>%
    flextable::set_table_properties(layout = "autofit") %>%
    flextable::set_header_labels(
        name = "Group", result = "IRITD3 Slope (95%CI)",
        "p_adjusted" = "FDR"
    ) %>%
    flextable::bold(i = NULL, part = "header") %>%
    flextable::fontsize(size = size_font, part = "all")

Group	IRITD3 Slope (95%CI)	FDR
Placebo	0.06 (-0.06 to 0.18)	0.3453
Felzartamab	-0.14 (-0.26 to -0.02)	0.0408
Intergroup Difference	-0.20 (-0.37 to -0.02)	0.0408

IRITD5

model_IRITD5 %>%
    dplyr::select(-p_value) %>%
    flextable::qflextable() %>%
    flextable::set_table_properties(layout = "autofit") %>%
    flextable::set_header_labels(
        name = "Group", result = "IRITD5 Slope (95%CI)",
        "p_adjusted" = "FDR"
    ) %>%
    flextable::bold(i = NULL, part = "header") %>%
    flextable::fontsize(size = size_font, part = "all")

Group	IRITD5 Slope (95%CI)	FDR
Placebo	0.09 (-0.04 to 0.22)	0.15720
Felzartamab	-0.13 (-0.26 to -0.01)	0.05895
Intergroup Difference	-0.23 (-0.41 to -0.05)	0.04200

2.5 PLOT THE RESULTS

# DEFINE LEGEND FOR VIOLIN PLOTS ####
panel_legend_violin <- plots_artANOVA %>%
    dplyr::filter(variable == "AMAT1") %>%
    pull(plot_violin) %>%
    ggpubr::get_legend() %>%
    ggpubr::as_ggplot() +
    theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))


# DEFINE LEGEND FOR REGRESSION PLOTS ####
panel_legend_regression <- plots_regression %>%
    dplyr::filter(variable == "IRRAT30") %>%
    pull(plot) %>%
    ggpubr::get_legend() %>%
    ggpubr::as_ggplot()


# DEFIN JOINT LEGEND FOR REGRESSION AND VIOLIN PLOTS ####
panel_legend_all_injury <- panel_legend_violin + panel_legend_regression


# DEFINE INJURY VIOLIN PLOTS ####
plots_violin <- plots_artANOVA %>%
    dplyr::filter(category %in% c("injury")) %>%
    pull(plot_violin)


# MAKE PANEL OF SLOPE AND VIOLIN PLOTS ####
plot_panels_all_injury <- patchwork::wrap_plots(
    plots_violin[[1]], plots_violin[[2]], plots_violin[[3]],
    plots_regression$plot[[1]], plots_regression$plot[[2]], plots_regression$plot[[3]],
    plots_regression$grob_table[[1]], plots_regression$grob_table[[2]], plots_regression$grob_table[[3]],
    nrow = 3,
    ncol = 3,
    guides = "collect",
    heights = c(1, 1, 0.5)
) +
    patchwork::plot_annotation(tag_levels = list(c(LETTERS[1:6], rep("", 3)))) &
    theme(
        legend.position = "none",
        # plot.title = element_blank(),
        axis.text = element_text(size = 10, colour = "black"),
        axis.title = element_text(size = 12, colour = "black"),
        plot.tag = element_text(size = 20, face = "bold", vjust = 1)
    )

ggarrange(
    panel_legend_all_injury,
    plot_panels_all_injury,
    nrow = 2,
    heights = c(0.15, 1.5)
)

3 Analysing the effect of treatment on gene expression

3.1 INTRODUCTION

This analysis assesses how felzartamab therapy affects gene expression measured in kidney allograft biopsies taken from patients with antibody mediated rejection. We use the Linear Models for Microarray Data (limma) pipeline for applying empirical Bayes mediated t-tests on individual genes. Specific contrasts to assess the interaction between treatment and time are specified to estimate the effect of treatment on each gene.

References:

• Ritchie, M.E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015).

3.2 DATA WRANGLING

First we need to load the necessary R packages and data.

# HOUSEKEEPING ####
# CRAN libraries
library(tidyverse) 
library(flextable) 
library(officer)
# Bioconductor libraries
library(Biobase) 
library(limma) 

## 
## Attaching package: 'limma'

## The following object is masked from 'package:BiocGenerics':
## 
##     plotMA

library(genefilter) 

## 
## Attaching package: 'genefilter'

## The following object is masked from 'package:patchwork':
## 
##     area

## The following object is masked from 'package:MASS':
## 
##     area

## The following object is masked from 'package:rstatix':
## 
##     Anova

## The following object is masked from 'package:readr':
## 
##     spec

# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/data/affymap219.RData")

# Seed for reproducibility
set.seed(42)

# DEFINE THE SET ####
set00 <- set[, set$Patient %nin% c(15, 18)]

• We filter the gene expression data by to remove genes with minimal variance across the population (IQR filtering).

# IQR FILTER THE DATA ####
f1 <- function(x) (IQR(x) > 0.5)
ff <- filterfun(f1)
if (!exists("selected")) {
    selected <- genefilter(set00, ff)
}
set01 <- set00[selected, ]

• We then retain a single probeset for each gene based on whichever has the highest mean expression in the population.

# KEEP UNIQUE GENES (keep probe with highest mean expression) ####
mean_exprs_by_probe <- set01 %>%
    exprs() %>%
    as_tibble(rownames = "AffyID") %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
    mutate(mean_exprs = set01 %>%
        exprs() %>% rowMeans(), .after = Symb)

genes <- mean_exprs_by_probe %>%
    group_by(Symb) %>%
    dplyr::slice_max(mean_exprs) %>%
    dplyr::filter(Symb != "") %>%
    distinct(Symb, .keep_all = TRUE) %>%
    pull(AffyID)

set <- set01[featureNames(set01) %in% genes, ]

3.3 REMOVE GENES THAT DIFFER AT BASELINE BETWEEN PLACEBO AND FELZARTAMAB ARMS

The felzartamab trial was a phase II randomized control trial (RCT). While RCTs are the ‘gold standard’ for accounting for confounders while establishing causality (i.e., the effect of treatment), they are not perfect, particularly when the number of study participants is low. Thus, we pre-screened gene expression data to remove any genes that differed between placebo and felzartamab arms at baseline. This was done to reduce the potential for interpreting the effect of treatment on genes with dissimilar baseline profiles.

# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set$Felzartamab_Followup %>% droplevels()


# BLOCK DESIGN week24 - baseline ####
design_block01 <- model.matrix(~ 0 + Felzartamab_Followup)
contrast_block_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupBaseline_Felzartamab) -(Felzartamab_FollowupBaseline_Placebo)",
    levels = design_block01
)

# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_block_1 <- limma::lmFit(set, design_block01)
cfit_block_1 <- limma::contrasts.fit(fit_block_1, contrast_block_01)
ebayes_block_1 <- limma::eBayes(cfit_block_1)
tab_block_1 <- limma::topTable(ebayes_block_1, adjust = "fdr", sort.by = "p", number = "all")

# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means_baseline <- fit_block_1 %>%
    avearrays() %>%
    data.frame() %>%
    rownames_to_column("AffyID") %>%
    tibble() %>%
    mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
    rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
    dplyr::select(AffyID, contains("Baseline"))

The top 20 genes that differ at baseline (sorted by p-value) between placebo and felzartamab arms are summarized below.

(Note that all FDR are > 0.05. This is a consequence of differential expression across a large selection of genes (5,267 in this case) with small samples sizes (i.e., n = 10 by treatment group). Nonetheless, genes with p < 0.05 will be excluded from subsequent analyses).

# FORMAT TOPTABLES ####
baseline_deg <- tab_block_1 %>%
    as_tibble(rownames = "AffyID") %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID") %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means_baseline, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb,
        Baseline_Placebo, Baseline_Felzartamab,
        t, logFC, P.Value, adj.P.Val,
    ) %>%
    dplyr::rename(p = P.Value, FDR = adj.P.Val)

baseline_deg %>%
    dplyr::slice_min(p, n = 20) %>%
    mutate_at(vars(t, logFC, FDR), ~ formatC(.x, format = "f", digits = 2)) %>%
    mutate_at(vars(p), ~ formatC(.x, format = "f", digits = 5)) %>%
    flextable::flextable()

AffyID	Symb	Baseline_Placebo	Baseline_Felzartamab	t	logFC	p	FDR
11727331_at	SPAG4	90	50	-4.23	-0.85	0.00008	0.33
11763725_a_at	HNRNPU	366	665	4.10	0.86	0.00013	0.33
11718286_a_at	UBL3	74	142	3.66	0.94	0.00053	0.45
11732590_x_at	TRIM5	42	63	3.65	0.59	0.00055	0.45
11725671_a_at	IPO7	42	66	3.54	0.64	0.00079	0.45
11741431_a_at	KCNJ16	203	311	3.49	0.61	0.00091	0.45
11726915_at	SLCO4C1	87	176	3.48	1.02	0.00096	0.45
11723109_a_at	GPBP1	153	263	3.46	0.78	0.00099	0.45
11739744_x_at	MXI1	123	191	3.45	0.64	0.00102	0.45
11720890_a_at	AGL	39	63	3.40	0.70	0.00121	0.45
11746431_a_at	FBXO11	28	45	3.38	0.68	0.00127	0.45
11751841_a_at	TNFRSF17	56	20	-3.36	-1.51	0.00134	0.45
11744367_a_at	ZMYND11	30	49	3.33	0.68	0.00150	0.45
11748857_a_at	AK3	371	578	3.32	0.64	0.00155	0.45
11727307_x_at	ZNF320	46	66	3.31	0.51	0.00160	0.45
11747648_a_at	TCAIM	96	152	3.29	0.66	0.00168	0.45
11728966_s_at	ZNF28	29	67	3.27	1.21	0.00181	0.45
11723387_a_at	PTAR1	149	261	3.25	0.81	0.00187	0.45
11739771_a_at	OXR1	31	51	3.24	0.73	0.00194	0.45
11730038_at	GNG13	44	28	-3.24	-0.63	0.00196	0.45

3.4 ESTIMATING TREATMENT EFFECT ON GENE EXPRESSION

We now carry out differential expression analysis to estimate the treatment effect on individual genes.

# DEFINE GENES SIMILAR AT BASELINE ####
genes_baseline <- baseline_deg %>%
    dplyr::filter(p > 0.05) %>%
    pull(AffyID)

# WRANGLE THE MEAN EXPRESSION DATA ####
mean_exprs_by_probe <- set01 %>%
    exprs() %>%
    as_tibble(rownames = "AffyID") %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
    mutate(
        mean_exprs = set01 %>%
            exprs() %>%
            rowMeans(),
        .after = Symb
    )

genes <- mean_exprs_by_probe %>%
    group_by(Symb) %>%
    dplyr::slice_max(mean_exprs) %>%
    dplyr::filter(Symb != "", AffyID %in% genes_baseline) %>%
    distinct(Symb, .keep_all = TRUE) %>%
    pull(AffyID)

set02 <- set01[featureNames(set01) %in% genes, ]

We need to define the model structure for the contrasts we intend to make. Here we will define 3 major contrast groups; 1) effect of time in placebo group, 2) effect of time in felzartamab group, and 3) the interactive effect of time and treatment. We do this for each time window (baseline to week 24, week 24 to week 52, and baseline to week 52). We also include a molecular estimate for % cortex as a covariate to adjust for difference in biopsy composition.

# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set02$Felzartamab_Followup %>% droplevels()
cortex <- set02$Cortexprob


# DESIGN ####
design <- model.matrix(~ 0 + Felzartamab_Followup + cortex)


# CONTRAST DESIGN week24 - baseline ####
contrast_interaction_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 -(Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)
contrast_felzartamab_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
    levels = design
)
contrast_placebo_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)

# CONTRAST DESIGN week52 - week24 ####
contrast_interaction_02 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
    levels = design
)
contrast_felzartamab_02 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2",
    levels = design
)
contrast_placebo_02 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
    levels = design
)


# CONTRAST DESIGN week52 - baseline####
contrast_interaction_03 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)
contrast_felzartamab_03 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
    levels = design
)
contrast_placebo_03 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)

We can check the specified contrast matrix for interaction terms. This is for the baseline to week 24 window.

contrast_interaction_01 %>% as_tibble(rownames = "level")

## # A tibble: 7 x 2
##   level                                    x =  (Felzartamab_FollowupWeek~1
##   <chr>                                                               <dbl>
## 1 Felzartamab_FollowupBaseline_Placebo                                  0.5
## 2 Felzartamab_FollowupWeek24_Placebo                                   -0.5
## 3 Felzartamab_FollowupWeek52_Placebo                                    0  
## 4 Felzartamab_FollowupBaseline_Felzartamab                             -0.5
## 5 Felzartamab_FollowupWeek24_Felzartamab                                0.5
## 6 Felzartamab_FollowupWeek52_Felzartamab                                0  
## 7 cortex                                                                0  
## # i abbreviated name:
## #   1: `x =  (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 -(Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2`

and this is for the baseline to week 52 window.

contrast_interaction_03 %>% as_tibble(rownames = "level")

## # A tibble: 7 x 2
##   level                                    x =  (Felzartamab_FollowupWeek~1
##   <chr>                                                               <dbl>
## 1 Felzartamab_FollowupBaseline_Placebo                                  0.5
## 2 Felzartamab_FollowupWeek24_Placebo                                    0  
## 3 Felzartamab_FollowupWeek52_Placebo                                   -0.5
## 4 Felzartamab_FollowupBaseline_Felzartamab                             -0.5
## 5 Felzartamab_FollowupWeek24_Felzartamab                                0  
## 6 Felzartamab_FollowupWeek52_Felzartamab                                0.5
## 7 cortex                                                                0  
## # i abbreviated name:
## #   1: `x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2`

We now fit the models we’ve specified

# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_interaction_1 <- limma::lmFit(set02, design)
cfit_interaction_1 <- limma::contrasts.fit(fit_interaction_1, contrast_interaction_01)
ebayes_interaction_1 <- limma::eBayes(cfit_interaction_1)
tab_interaction_1 <- limma::topTable(ebayes_interaction_1, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_interaction_1$stdev.unscaled * sqrt(ebayes_interaction_1$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_felzartamab_1 <- limma::lmFit(set02, design)
cfit_felzartamab_1 <- limma::contrasts.fit(fit_felzartamab_1, contrast_felzartamab_01)
ebayes_felzartamab_1 <- limma::eBayes(cfit_felzartamab_1)
tab_felzartamab_1 <- limma::topTable(ebayes_felzartamab_1, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_felzartamab_1$stdev.unscaled * sqrt(ebayes_felzartamab_1$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_placebo_1 <- limma::lmFit(set02, design)
cfit_placebo_1 <- limma::contrasts.fit(fit_placebo_1, contrast_placebo_01)
ebayes_placebo_1 <- limma::eBayes(cfit_placebo_1)
tab_placebo_1 <- limma::topTable(ebayes_placebo_1, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_placebo_1$stdev.unscaled * sqrt(ebayes_placebo_1$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")


# FIT BLOCK week52 - week24 LIMMA MODEL ####
fit_interaction_2 <- limma::lmFit(set02, design)
cfit_interaction_2 <- limma::contrasts.fit(fit_interaction_2, contrast_interaction_02)
ebayes_interaction_2 <- limma::eBayes(cfit_interaction_2)
tab_interaction_2 <- limma::topTable(ebayes_interaction_2, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_interaction_2$stdev.unscaled * sqrt(ebayes_interaction_2$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_felzartamab_2 <- limma::lmFit(set02, design)
cfit_felzartamab_2 <- limma::contrasts.fit(fit_felzartamab_2, contrast_felzartamab_02)
ebayes_felzartamab_2 <- limma::eBayes(cfit_felzartamab_2)
tab_felzartamab_2 <- limma::topTable(ebayes_felzartamab_2, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_felzartamab_2$stdev.unscaled * sqrt(ebayes_felzartamab_2$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_placebo_2 <- limma::lmFit(set02, design)
cfit_placebo_2 <- limma::contrasts.fit(fit_placebo_2, contrast_placebo_02)
ebayes_placebo_2 <- limma::eBayes(cfit_placebo_2)
tab_placebo_2 <- limma::topTable(ebayes_placebo_2, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_placebo_2$stdev.unscaled * sqrt(ebayes_placebo_2$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")


# FIT BLOCK week52 - baseline LIMMA MODEL ####
fit_interaction_3 <- limma::lmFit(set02, design)
cfit_interaction_3 <- limma::contrasts.fit(fit_interaction_3, contrast_interaction_03)
ebayes_interaction_3 <- limma::eBayes(cfit_interaction_3)
tab_interaction_3 <- limma::topTable(ebayes_interaction_3, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_interaction_3$stdev.unscaled * sqrt(ebayes_interaction_3$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_felzartamab_3 <- limma::lmFit(set02, design)
cfit_felzartamab_3 <- limma::contrasts.fit(fit_felzartamab_3, contrast_felzartamab_03)
ebayes_felzartamab_3 <- limma::eBayes(cfit_felzartamab_3)
tab_felzartamab_3 <- limma::topTable(ebayes_felzartamab_3, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_felzartamab_3$stdev.unscaled * sqrt(ebayes_felzartamab_3$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_placebo_3 <- limma::lmFit(set02, design)
cfit_placebo_3 <- limma::contrasts.fit(fit_placebo_3, contrast_placebo_03)
ebayes_placebo_3 <- limma::eBayes(cfit_placebo_3)
tab_placebo_3 <- limma::topTable(ebayes_placebo_3, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_placebo_3$stdev.unscaled * sqrt(ebayes_placebo_3$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

3.5 ORGANIZE THE RESULTS

The results are now organized for tabulation

# MERGE BLOCKS ####
tab_block_1 <- tab_interaction_1 %>%
    left_join(
        tab_felzartamab_1 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(flogFC = logFC),
        by = "AffyID"
    ) %>%
    left_join(
        tab_placebo_1 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(plogFC = logFC),
        by = "AffyID"
    ) %>%
    relocate(plogFC, flogFC, .before = logFC) %>%
    arrange(P.Value)

tab_block_2 <- tab_interaction_2 %>%
    left_join(
        tab_felzartamab_2 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(flogFC = logFC),
        by = "AffyID"
    ) %>%
    left_join(
        tab_placebo_2 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(plogFC = logFC),
        by = "AffyID"
    ) %>%
    relocate(plogFC, flogFC, .before = logFC) %>%
    arrange(P.Value)

tab_block_3 <- tab_interaction_3 %>%
    left_join(
        tab_felzartamab_3 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(flogFC = logFC),
        by = "AffyID"
    ) %>%
    left_join(
        tab_placebo_3 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(plogFC = logFC),
        by = "AffyID"
    ) %>%
    relocate(plogFC, flogFC, .before = logFC) %>%
    arrange(P.Value)


# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means <- fit_interaction_1 %>%
    avearrays() %>%
    as_tibble(rownames = "AffyID") %>%
    mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
    rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
    dplyr::select(-contains("Week12"), -any_of(c("cortex")))


# FORMAT TOPTABLES ####
table_interaction_1 <- tab_block_1 %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb, Gene, PBT,
        t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
        all_of(colnames(means))
    ) %>%
    dplyr::rename(
        p = P.Value,
        FDR = adj.P.Val
    )

table_interaction_2 <- tab_block_2 %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb, Gene, PBT,
        t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
        all_of(colnames(means))
    ) %>%
    dplyr::rename(
        p = P.Value,
        FDR = adj.P.Val
    )

table_interaction_3 <- tab_block_3 %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb, Gene, PBT,
        t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
        all_of(colnames(means))
    ) %>%
    dplyr::rename(
        p = P.Value,
        FDR = adj.P.Val
    )

deg <- tibble(
    design = c(
        "Baseline_vs_Week24",
        "Week24_vs_Week52",
        "Baseline_vs_Week52"
    ),
    table = list(
        table_interaction_1,
        table_interaction_2,
        table_interaction_3
    )
)

3.6 SUMMARIZE THE RESULTS

# FORMAT TABLES TO MAKE FLEXTABLES ####
table_deg <- deg %>%
    expand_grid(direction = c("all", "increased", "decreased")) %>%
    relocate(direction, .after = design) %>%
    mutate(
        gene_tables = pmap(
            list(direction, table),
            function(direction, table) {
                if (direction == "increased") {
                    table <- table %>%
                        dplyr::filter(logFC > 0) %>%
                        arrange(p)
                } else if (direction == "decreased") {
                    table <- table %>%
                        dplyr::filter(logFC < 0) %>%
                        arrange(p)
                }
                table %>%
                    dplyr::select(-t, -se, -FDR, -contains("AffyID"), -contains("MMDx")) %>%
                    dplyr::slice(1:20) %>%
                    mutate(
                        Gene = Gene %>% str_remove("///.*"),
                        plogFC = plogFC %>% round(2),
                        flogFC = flogFC %>% round(2),
                        logFC = logFC %>% round(2),
                        p = case_when(
                            p < 0.0001 ~ p %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p %>% formatC(digits = 4, format = "f")
                        )
                    )
            }
        )
    )


# GLOBAL PARAMETERS FOR FLEXTABLES ####
header1 <- c(
    "Gene\nsymbol", "Gene", "PBT",
    rep("Differential expression", 4),
    rep("Mean expression by group", 6)
)
header2 <- c(
    "Gene\nsymbol", "Gene", "PBT",
    "\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
    "\u394\u394\nlogFC", "\u394\u394\nP",
    rep("Placebo", 3), rep("Felzartamab", 3)
)
header3 <- c(
    "Gene\nsymbol", "Gene", "PBT",
    "\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
    "\u394\u394\nlogFC", "\u394\u394\nP",
    "Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)", "Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)"
)

# MAKE FORMATTED FLEXTABLES ####
flextables <- table_deg %>%
    mutate(
        flextable = pmap(
            list(design, gene_tables),
            function(design, gene_tables) {
                if (design == "Baseline_vs_Week24") {
                    title <- paste("Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
                } else if (design == "Week24_vs_Week52") {
                    title <- paste("Table i. Top 20 differentially expressed genes between week24 and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
                } else if (design == "Baseline_vs_Week52") {
                    title <- paste("Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
                }
                gene_tables %>%
                    arrange(logFC) %>%
                    flextable::flextable() %>%
                    flextable::delete_part("header") %>%
                    flextable::add_header_row(top = TRUE, values = header3) %>%
                    flextable::add_header_row(top = TRUE, values = header2) %>%
                    flextable::add_header_row(top = TRUE, values = header1) %>%
                    flextable::add_header_row(top = TRUE, values = rep(title, ncol_keys(.))) %>%
                    flextable::merge_v(part = "header") %>%
                    flextable::merge_h(part = "header") %>%
                    flextable::border_remove() %>%
                    flextable::border(part = "header", border = officer::fp_border()) %>%
                    flextable::border(part = "body", border = officer::fp_border()) %>%
                    flextable::align(align = "center") %>%
                    flextable::align(align = "center", part = "header") %>%
                    flextable::font(fontname = "Arial", part = "all") %>%
                    flextable::fontsize(size = 7, part = "all") %>%
                    flextable::fontsize(size = 7, part = "footer") %>%
                    flextable::fontsize(i = 1, size = 12, part = "header") %>%
                    flextable::bold(part = "header") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::padding(padding = 0, part = "all") %>%
                    flextable::autofit()
            }
        )
    )

Genes suppressed by felzartamab in week 24 compared to baseline
(Note that all ΔΔFDR are > 0.05. This is a consequence of differential expression across a large selection of genes (5,267 in this case) with small samples sizes (i.e., n = 10 by treatment group). Thus, treatment effects on individual genes should not be interpreted. Instead, geneset enrichment analysis (GSEA) will be used to assess bulk trends).

flextables %>%
    dplyr::filter(design == "Baseline_vs_Week24", direction == "decreased") %>%
    pull(flextable) %>%
    pluck(1)

Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value)
Gene symbol	Gene	PBT	Differential expression				Mean expression by group
			Δ placebo logFC	Δ felzartamab logFC	ΔΔ logFC	ΔΔ P	Placebo			Felzartamab
							Baseline (N=10)	Week24 (N=10)	Week52 (N=10)	Baseline (N=10)	Week24 (N=10)	Week52 (N=10)
CXCL11	chemokine (C-X-C motif) ligand 11	ABMR-RAT,GRIT3,RAT,Rej-RAT	-0.16	-1.01	-0.85	0.0176	494	395	372	617	153	317
CXCL9	chemokine (C-X-C motif) ligand 9	ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT	-0.14	-0.94	-0.79	0.0290	1,344	1,102	1,384	1,684	459	1,036
CXCL10	chemokine (C-X-C motif) ligand 10	ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT	-0.21	-0.98	-0.77	0.0228	844	632	720	1,063	274	615
IDO1	indoleamine 2,3-dioxygenase 1	ABMR-RAT,GRIT3,RAT,Rej-RAT	-0.09	-0.77	-0.68	0.0103	338	298	296	409	140	251
LYZ	lysozyme	IRITD5,QCMAT	-0.01	-0.67	-0.66	0.0386	1,676	1,658	1,884	2,244	884	1,184
FCGR3A	Fc fragment of IgG, low affinity IIIa, receptor (CD16a)	ABMR-RAT,GRIT2,IRRAT950,RAT,Rej-RAT	0.06	-0.58	-0.64	0.0167	550	602	613	571	256	352
GZMB	granzyme B	ABMR-RAT,QCAT,RAT,Rej-RAT	-0.01	-0.62	-0.61	0.0019	149	147	130	143	61	114
GBP1	guanylate binding protein 1, interferon-inducible	ABMR-RAT,GRIT3,RAT,Rej-RAT	-0.03	-0.61	-0.59	0.0096	678	655	643	909	389	556
GNLY	granulysin	ABMR-RAT,DSAST,QCAT,RAT,Rej-RAT	-0.18	-0.75	-0.57	0.0285	192	149	134	165	58	100
SH2D1B	SH2 domain containing 1B	ABMR-RAT,DSAST,NKB,RAT	0.07	-0.44	-0.51	0.0010	38	42	34	40	22	33
FCGR1A	Fc fragment of IgG, high affinity Ia, receptor (CD64)	GRIT2,IRRAT950,RAT,TCMR-RAT	0.09	-0.41	-0.50	0.0304	131	148	161	152	86	86
KLRD1	killer cell lectin-like receptor subfamily D, member 1	ABMR-RAT,RAT,Rej-RAT	-0.06	-0.49	-0.43	0.0295	110	101	96	126	63	102
EPB41L3	erythrocyte membrane protein band 4.1-like 3		0.22	-0.20	-0.42	0.0388	209	286	258	339	258	283
LOC339059	uncharacterized LOC339059		0.26	-0.15	-0.41	0.0203	51	73	55	66	53	59
APOL1	apolipoprotein L1	ABMR-RAT,GRIT3,RAT,Rej-RAT	0.07	-0.32	-0.39	0.0205	626	687	744	754	483	636
LYPD5	LY6/PLAUR domain containing 5	ABMR-RAT,RAT,Rej-RAT	0.04	-0.33	-0.37	0.0255	27	28	27	30	19	26
TRAV12-2	T cell receptor alpha variable 12-2		0.25	-0.10	-0.35	0.0173	39	55	39	49	43	43
IL18RAP	interleukin 18 receptor accessory protein		0.02	-0.32	-0.34	0.0233	38	39	41	42	27	33
PRF1	perforin 1 (pore forming protein)	ABMR-RAT,QCAT,RAT,Rej-RAT	-0.02	-0.35	-0.33	0.0382	152	148	137	145	89	123
PMM2	phosphomannomutase 2		0.20	-0.07	-0.27	0.0380	38	50	41	46	42	49

Genes suppressed by felzartamab in week 52 compared to baseline
(Note that all ΔΔFDR are > 0.05. This is a consequence of differential expression across a large selection of genes (5,267 in this case) with small samples sizes (i.e., n = 10 by treatment group). Thus, treatment effects on individual genes should not be interpreted. Instead, geneset enrichment analysis (GSEA) will be used to assess bulk trends).

flextables %>%
    dplyr::filter(design == "Baseline_vs_Week52", direction == "decreased") %>%
    pull(flextable) %>%
    pluck(1)

Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)
Gene symbol	Gene	PBT	Differential expression				Mean expression by group
			Δ placebo logFC	Δ felzartamab logFC	ΔΔ logFC	ΔΔ P	Placebo			Felzartamab
							Baseline (N=10)	Week24 (N=10)	Week52 (N=10)	Baseline (N=10)	Week24 (N=10)	Week52 (N=10)
HNRNPK	heterogeneous nuclear ribonucleoprotein K		0.28	-0.38	-0.66	0.0055	818	1,044	1,203	1,096	1,059	649
COL1A1	collagen, type I, alpha 1	COLA1356,FICOL,IRITD5,LivGST_UP	0.38	-0.26	-0.64	0.0075	310	333	526	329	320	230
NFIX	nuclear factor I/X (CCAAT-binding transcription factor)	CT1	0.42	-0.15	-0.57	0.0012	162	179	289	221	167	180
RHOA	ras homolog family member A	IRRAT950	0.28	-0.29	-0.57	0.0049	691	803	1,017	923	809	620
COL1A2	collagen, type I, alpha 2	FICOL,IRITD5,LivGST_UP	0.26	-0.31	-0.57	0.0118	1,238	1,152	1,769	1,534	1,194	997
ITGB3	integrin beta 3	IRRAT30,IRRAT950	0.34	-0.21	-0.55	0.0100	105	137	169	110	103	82
LOC100129518	uncharacterized LOC100129518	GRIT3,IRRAT950	0.25	-0.22	-0.47	0.0050	429	546	605	527	461	390
ACTB	actin, beta	IRRAT950	0.15	-0.31	-0.47	0.0077	6,451	7,352	7,993	7,274	5,662	4,720
SPARC	secreted protein, acidic, cysteine-rich (osteonectin)	IRITD3	0.16	-0.30	-0.46	0.0100	2,411	2,532	3,011	2,343	2,094	1,550
DNAJA4	DnaJ (Hsp40) homolog, subfamily A, member 4	IRRAT950	0.27	-0.14	-0.42	0.0095	49	68	72	68	61	56
SLFN5	schlafen family member 5		0.19	-0.20	-0.40	0.0050	39	44	51	49	38	37
CYR61	cysteine-rich, angiogenic inducer, 61	IRITD3,IRRAT950	0.14	-0.26	-0.40	0.0074	271	268	330	338	234	236
SDCBP	syndecan binding protein		0.24	-0.14	-0.38	0.0120	1,135	1,404	1,582	1,478	1,358	1,224
SBSPON	somatomedin B and thrombospondin type 1 domain containing		0.31	-0.08	-0.38	0.0123	217	250	332	188	249	169
PAFAH1B2	platelet-activating factor acetylhydrolase 1b, catalytic subunit 2 (30kDa)		0.23	-0.13	-0.36	0.0099	309	372	423	386	396	323
TMEM2	transmembrane protein 2		0.10	-0.26	-0.36	0.0118	270	307	311	318	259	222
CALM2	calmodulin 2 (phosphorylase kinase, delta)	CT1,IRITD3	0.12	-0.25	-0.36	0.0138	1,304	1,427	1,538	1,623	1,342	1,155
WNK1	WNK lysine deficient protein kinase 1		0.26	-0.08	-0.34	0.0104	703	884	1,012	894	996	799
CEBPB	CCAAT/enhancer binding protein (C/EBP), beta	IRITD3	0.11	-0.22	-0.33	0.0123	419	452	488	417	333	307
STAT6	signal transducer and activator of transcription 6, interleukin-4 induced		0.09	-0.20	-0.30	0.0144	671	710	765	746	616	563

4 Geneset enrichment analysis to identify pathways affected by felzartamab treatment

4.1 INTRODUCTION

Geneset enrichment analysis (GSEA) is applied to identify pathways affected by felzartamab treatment. Multiple libraries were screened to encompass a broad spectrum of biological pathways, including Gene Ontology (GO), Disease Ontology Semantic and Enrichment analysis (DOSE), Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome, Molecular Signatures Database (MsigDB), and WikiPathways (wiki) libraries. An additional de-novo injury geneset was also assesed. This injury geneset was based on single-nuclei transcriptomics assessed from healthy and injury kidney allografts. Geneset enrichment was chosen specifcally because it is a rank-based test and does not require cutoffs based on p-values or t-values. Nonetheless, here we decided apriori that we would pre-screen genes based on their p-values from differential expression analyses described in Section 3. Note that there are many ways for carrying out pathway and functional enrichment analyses. Also note that gene libraries are frequently updated and you may find slightly different results based on the packge version being used. For example, here we encountered such a difference in the DOSE library, where versions later than 3.30.1 identified two different pathways. Thus, we have included a tarball of the package version used for the original analyses in this paper for convenient installation and replication of the presented findings.

References:

• Yu, G., Wang, L.G., Han, Y. & He, Q.Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284-287 (2012).

• Wu, T., et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2, 100141 (2021).

• Hinze, C., et al. Single-cell transcriptomics reveals common epithelial response patterns in human acute kidney injury. Genome Med 14, 103 (2022).

4.2 DATA WRANGLING

First we need to load the necessary R packages and data.

# HOUSEKEEPING ####
# CRAN libraries
library(tidyverse)
library(flextable) 
library(officer) 
library(msigdbr)
# Bioconductor libraries
# DOSE v 3.30.1 required for reproducibility
# install.packages("natmed/data/DOSE_3.30.1.tar.gz", repos = NULL, type = "source")
library(DOSE)
library(ReactomePA)
library(Biobase) 
library(clusterProfiler)
library(org.Hs.eg.db) 
library(pRoloc) 

## Warning in fun(libname, pkgname): mzR has been built against a different Rcpp version (1.0.12)
## than is installed on your system (1.0.13). This might lead to errors
## when loading mzR. If you encounter such issues, please send a report,
## including the output of sessionInfo() to the Bioc support forum at 
## https://support.bioconductor.org/. For details see also
## https://github.com/sneumann/mzR/wiki/mzR-Rcpp-compiler-linker-issue.

# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# load DE data
load("C:/R/CD38-effect-of-treatment/natmed/results/deg.RData")

# Seed for reproducibility
set.seed(42)

4.3 GENESET ENRICHMENT ANALYSES

Geneset enrichment analyses using the ClusterProfiler package have simlar, but not identical implementation. This section will carry out GSEA on all described libraries, after which the results will be collated and summarized.

GO library

# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
    mutate(
        genes_gsea = map(
            table,
            function(table) {
                table %>%
                    dplyr::filter(p < 0.05) %>%
                    arrange(t %>% dplyr::desc()) %>%
                    dplyr::mutate(
                        ENTREZID = bitr(Symb,
                            fromType = "SYMBOL",
                            toType = c("ENTREZID"),
                            OrgDb = org.Hs.eg.db,
                            drop = FALSE
                        ) %>% pull(ENTREZID)
                    ) %>%
                    drop_na(ENTREZID) %>%
                    dplyr::select(t, ENTREZID) %>%
                    pull(t, ENTREZID)
            }
        )
    )


# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_go <- data %>%
    mutate(
        gsea = map(
            genes_gsea,
            function(genes_gsea) {
                clusterProfiler::gseGO(
                    gene = genes_gsea, ont = "BP", OrgDb = org.Hs.eg.db,
                    minGSSize = 10, maxGSSize = 200,
                    pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
                ) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
            }
        )
    )


# SIMPLIFY GO TERMS ####
set.seed(42)
gsea_go_simplified <- gsea_go %>%
    mutate(gsea_simplified = map(gsea, clusterProfiler::simplify))


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_go_simplified_formatted <- gsea_go_simplified %>%
    mutate(
        gsea_tables = map(
            gsea_simplified,
            function(gsea_simplified) {
                gsea_simplified %>%
                    as_tibble() %>%
                    arrange(pvalue) %>%
                    mutate(
                        sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
                        Description = factor(Description, levels = Description, ordered = TRUE)
                    ) 
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_go_tables <- gsea_go_simplified_formatted %>%
    mutate(
        gsea_flextables = map(
            gsea_tables,
            function(gsea_tables) {
                gsea_tables %>%
                    slice_min(pvalue, n = 20) %>%
                    mutate(
                        NES = NES %>% round(2),
                        pvalue = case_when(
                            pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
                            TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
                        ),
                        FDR = case_when(
                            p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
                        ),
                    ) %>%
                    dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
                    flextable::flextable() %>%
                    flextable::border_remove() %>%
                    flextable::border(border = officer::fp_border(), part = "all") %>%
                    flextable::align(align = "center", part = "all") %>%
                    flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::autofit()
            }
        )
    )

# PREPARE THE RESULTS FOR EXPORT ####
gsea_go <- gsea_go_tables %>%
    mutate(db = "GO", .before = 1)
names(gsea_go$gsea_flextables) <- gsea_go$design

DOSE library

# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
    mutate(
        genes_gsea = map(
            table,
            function(table) {
                table %>%
                    dplyr::filter(p < 0.05) %>%
                    arrange(t %>% dplyr::desc()) %>%
                    dplyr::mutate(
                        ENTREZID = clusterProfiler::bitr(Symb,
                            fromType = "SYMBOL",
                            toType = c("ENTREZID"),
                            OrgDb = org.Hs.eg.db,
                            drop = FALSE
                        ) %>% pull(ENTREZID)
                    ) %>%
                    drop_na(ENTREZID) %>%
                    dplyr::select(t, ENTREZID) %>%
                    pull(t, ENTREZID)
            }
        )
    )


# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_dose <- data %>%
    mutate(
        gsea = map(
            genes_gsea,
            function(genes_gsea) {
                DOSE::gseDO(
                    gene = genes_gsea,
                    minGSSize = 10, maxGSSize = 200,
                    pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
                ) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_dose_formatted <- gsea_dose %>%
    mutate(
        gsea_tables = map(
            gsea,
            function(gsea) {
                gsea %>%
                    as_tibble() %>%
                    arrange(pvalue) %>%
                    mutate(sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated")) 
            }
        )
    )



# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_dose_tables <- gsea_dose_formatted %>%
    mutate(
        gsea_flextables = map(
            gsea_tables,
            function(gsea_tables) {
                gsea_tables %>%
                    slice_min(pvalue, n = 20) %>%
                    mutate(
                        NES = NES %>% round(2),
                        pvalue = case_when(
                            pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
                            TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
                        ),
                        FDR = case_when(
                            p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
                        ),
                    ) %>%
                    dplyr::select(Description, setSize, NES, pvalue, FDR,  dplyr::any_of("core_enrichment")) %>%
                    flextable::flextable() %>%
                    flextable::border_remove() %>%
                    flextable::border(border = officer::fp_border(), part = "all") %>%
                    flextable::align(align = "center", part = "all") %>%
                    flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::autofit()
            }
        )
    )

# PREPARE THE RESULTS FOR EXPORT ####
gsea_dose <- gsea_dose_tables %>%
    mutate(db = "DOSE", .before = 1)
names(gsea_dose$gsea_flextables) <- gsea_dose$design

KEGG library

# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
    mutate(
        genes_gsea = map(
            table,
            function(table) {
                table %>%
                    dplyr::filter(p < 0.05) %>%
                    arrange(t %>% dplyr::desc()) %>%
                    dplyr::mutate(
                        ENTREZID = bitr(Symb,
                            fromType = "SYMBOL",
                            toType = c("ENTREZID"),
                            OrgDb = org.Hs.eg.db,
                            drop = FALSE
                        ) %>% pull(ENTREZID)
                    ) %>%
                    drop_na(ENTREZID) %>%
                    dplyr::select(t, ENTREZID) %>%
                    pull(t, ENTREZID)
            }
        )
    )


# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_kegg <- data %>%
    mutate(
        gsea = map(
            genes_gsea,
            function(genes_gsea) {
                clusterProfiler::gseKEGG(
                    gene = genes_gsea, organism = "hsa",
                    minGSSize = 10, maxGSSize = 200,
                    pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
                ) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_kegg_formatted <- gsea_kegg %>%
    mutate(
        gsea_tables = map(
            gsea,
            function(gsea) {
                gsea %>%
                    as_tibble() %>%
                    arrange(pvalue) %>%
                    mutate(
                        sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
                        Description = factor(Description, levels = Description, ordered = TRUE)
                    )
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_kegg_tables <- gsea_kegg_formatted %>%
    mutate(
        gsea_flextables = map(
            gsea_tables,
            function(gsea_tables) {
                gsea_tables %>%
                    slice_min(pvalue, n = 20) %>%
                    mutate(
                        NES = NES %>% round(2),
                        pvalue = case_when(
                            pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
                            TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
                        ),
                        FDR = case_when(
                            p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
                        ),
                    ) %>%
                    dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
                    flextable::flextable() %>%
                    flextable::border_remove() %>%
                    flextable::border(border = officer::fp_border(), part = "all") %>%
                    flextable::align(align = "center", part = "all") %>%
                    flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::autofit()
            }
        )
    )


# PREPARE THE RESULTS FOR EXPORT ####
gsea_kegg <- gsea_kegg_tables %>%
    mutate(db = "kegg", .before = 1)
names(gsea_kegg$gsea_flextables) <- gsea_kegg$design

MSigDB library

# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
    mutate(
        genes_gsea = map(
            table,
            function(table) {
                table %>%
                    dplyr::filter(p < 0.05) %>%
                    arrange(t %>% dplyr::desc()) %>%
                    dplyr::mutate(
                        ENTREZID = bitr(Symb,
                            fromType = "SYMBOL",
                            toType = c("ENTREZID"),
                            OrgDb = org.Hs.eg.db,
                            drop = FALSE
                        ) %>% pull(ENTREZID)
                    ) %>%
                    drop_na(ENTREZID) %>%
                    dplyr::select(t, ENTREZID) %>%
                    pull(t, ENTREZID)
            }
        )
    )


# DEFINE MSigDB PATHWAYS ####
gene_sets <- msigdbr(category = "H")
map <- gene_sets[, c("gs_name", "entrez_gene")]
map$entrez_gene <- as.character(map$entrez_gene)


# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_msigdb <- data %>%
    mutate(
        gsea = map(
            genes_gsea,
            function(genes_gsea) {
                clusterProfiler::GSEA(
                    gene = genes_gsea, TERM2GENE = map,
                    minGSSize = 10, maxGSSize = 200,
                    pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
                ) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_msigdb_formatted <- gsea_msigdb %>%
    mutate(
        gsea_tables = map(
            gsea,
            function(gsea) {
                gsea %>%
                    as_tibble() %>%
                    arrange(pvalue) %>%
                    mutate(
                        sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
                        Description = factor(Description, levels = Description, ordered = TRUE)
                    ) 
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_msigdb_tables <- gsea_msigdb_formatted %>%
    mutate(
        gsea_flextables = map(
            gsea_tables,
            function(gsea_tables) {
                gsea_tables %>%
                    # slice_min(pvalue, n = 20) %>%
                    mutate(
                        NES = NES %>% round(2),
                        pvalue = case_when(
                            pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
                            TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
                        ),
                        FDR = case_when(
                            p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
                        ),
                    ) %>%
                    dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
                    flextable::flextable() %>%
                    flextable::border_remove() %>%
                    flextable::border(border = officer::fp_border(), part = "all") %>%
                    flextable::align(align = "center", part = "all") %>%
                    flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::autofit()
            }
        )
    )


# PREPARE THE RESULTS FOR EXPORT ####
gsea_msigdb <- gsea_msigdb_tables %>%
    mutate(db = "msigdb", .before = 1)
names(gsea_msigdb$gsea_flextables) <- gsea_msigdb$design

REACTOME library

# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
    mutate(
        genes_gsea = map(
            table,
            function(table) {
                table %>%
                    dplyr::filter(p < 0.05) %>%
                    arrange(t %>% dplyr::desc()) %>%
                    dplyr::mutate(
                        ENTREZID = bitr(Symb,
                            fromType = "SYMBOL",
                            toType = c("ENTREZID"),
                            OrgDb = org.Hs.eg.db,
                            drop = FALSE
                        ) %>% pull(ENTREZID)
                    ) %>%
                    drop_na(ENTREZID) %>%
                    dplyr::select(t, ENTREZID) %>%
                    pull(t, ENTREZID)
            }
        )
    )


# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_reactome <- data %>%
    mutate(
        gsea = map(
            genes_gsea,
            function(genes_gsea) {
                ReactomePA::gsePathway(
                    gene = genes_gsea, organism = "human",
                    minGSSize = 10, maxGSSize = 200,
                    pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
                ) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_reactome_formatted <- gsea_reactome %>%
    mutate(
        gsea_tables = map(
            gsea,
            function(gsea) {
                gsea %>%
                    as_tibble() %>%
                    arrange(pvalue) %>%
                    mutate(
                        sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
                        Description = factor(Description, levels = Description, ordered = TRUE)
                    ) 
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_reactome_tables <- gsea_reactome_formatted %>%
    mutate(
        gsea_flextables = map(
            gsea_tables,
            function(gsea_tables) {
                gsea_tables %>%
                    slice_min(pvalue, n = 20) %>%
                    mutate(
                        NES = NES %>% round(2),
                        pvalue = case_when(
                            pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
                            TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
                        ),
                        FDR = case_when(
                            p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
                        ),
                    ) %>%
                    dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
                    flextable::flextable() %>%
                    flextable::border_remove() %>%
                    flextable::border(border = officer::fp_border(), part = "all") %>%
                    flextable::align(align = "center", part = "all") %>%
                    flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::autofit()
            }
        )
    )


# PREPARE THE RESULTS FOR EXPORT ####
gsea_reactome <- gsea_reactome_tables %>%
    mutate(db = "Reactome", .before = 1)
names(gsea_reactome$gsea_flextables) <- gsea_reactome$design

wikiPathways library

# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
    mutate(
        genes_gsea = map(
            table,
            function(table) {
                table %>%
                    dplyr::filter(p < 0.05) %>%
                    arrange(t %>% dplyr::desc()) %>%
                    dplyr::mutate(
                        ENTREZID = bitr(Symb,
                            fromType = "SYMBOL",
                            toType = c("ENTREZID"),
                            OrgDb = org.Hs.eg.db,
                            drop = FALSE
                        ) %>% pull(ENTREZID)
                    ) %>%
                    drop_na(ENTREZID) %>%
                    dplyr::select(t, ENTREZID) %>%
                    pull(t, ENTREZID)
            }
        )
    )


# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_wiki <- data %>%
    mutate(
        gsea = map(
            genes_gsea,
            function(genes_gsea) {
                clusterProfiler::gseWP(
                    gene = genes_gsea, organism = "Homo sapiens",
                    minGSSize = 10, maxGSSize = 200,
                    pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
                ) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_wiki_formatted <- gsea_wiki %>%
    mutate(
        gsea_tables = map(
            gsea,
            function(gsea) {
                gsea %>%
                    as_tibble() %>%
                    arrange(pvalue) %>%
                    mutate(
                        sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
                        Description = factor(Description, levels = Description, ordered = TRUE)
                    ) 
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_wiki_tables <- gsea_wiki_formatted %>%
    mutate(
        gsea_flextables = map(
            gsea_tables,
            function(gsea_tables) {
                gsea_tables %>%
                    # slice_min(pvalue, n = 20) %>%
                    mutate(
                        NES = NES %>% round(2),
                        pvalue = case_when(
                            pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
                            TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
                        ),
                        FDR = case_when(
                            p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
                        ),
                    ) %>%
                    dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
                    flextable::flextable() %>%
                    flextable::border_remove() %>%
                    flextable::border(border = officer::fp_border(), part = "all") %>%
                    flextable::align(align = "center", part = "all") %>%
                    flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::autofit()
            }
        )
    )


# PREPARE THE RESULTS FOR EXPORT ####
gsea_wiki <- gsea_wiki_tables %>%
    mutate(db = "wiki", .before = 1)
names(gsea_wiki$gsea_flextables) <- gsea_wiki$design

AKI geneset

# load gene lists
load("C:/R/CD38-effect-of-treatment/natmed/data/hinze2022_injury_markers.RData")


# DEFINE INJURY PATHWAYS ####
map_aki <- injury_markers %>% dplyr::select(gs_name, entrez_gene)
map_aki_joined <- injury_markers %>%
    dplyr::select(gs_name, entrez_gene) %>%
    mutate(gs_name = "AKI")


# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
    mutate(
        genes_gsea = map(
            table,
            function(table) {
                table %>%
                    dplyr::filter(p < 0.05) %>%
                    arrange(t %>% dplyr::desc()) %>%
                    dplyr::mutate(
                        ENTREZID = bitr(Symb,
                            fromType = "SYMBOL",
                            toType = c("ENTREZID"),
                            OrgDb = org.Hs.eg.db,
                            drop = FALSE
                        ) %>% pull(ENTREZID)
                    ) %>%
                    drop_na(ENTREZID) %>%
                    dplyr::select(t, ENTREZID) %>%
                    pull(t, ENTREZID)
            }
        )
    )


# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_aki <- data %>%
    mutate(
        gsea = map(
            genes_gsea,
            function(genes_gsea) {
                clusterProfiler::GSEA(
                    gene = genes_gsea, TERM2GENE = map_aki_joined,
                    minGSSize = 5, maxGSSize = Inf,
                    pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
                ) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_aki_formatted <- gsea_aki %>%
    mutate(
        gsea_tables = map(
            gsea,
            function(gsea) {
                gsea %>%
                    as_tibble() %>%
                    arrange(pvalue) %>%
                    mutate(
                        sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
                        Description = "Cellular response to acute kidney injury"
                    )
            }
        )
    )


# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_aki_tables <- gsea_aki_formatted %>%
    mutate(
        gsea_flextables = map(
            gsea_tables,
            function(gsea_tables) {
                gsea_tables %>%
                    mutate(
                        NES = NES %>% round(2),
                        pvalue = case_when(
                            pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
                            TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
                        ),
                        FDR = case_when(
                            p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
                        ),
                    ) %>%
                    dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
                    flextable::flextable() %>%
                    flextable::border_remove() %>%
                    flextable::border(border = officer::fp_border(), part = "all") %>%
                    flextable::align(align = "center", part = "all") %>%
                    flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::autofit()
            }
        )
    )


# PREPARE THE RESULTS FOR EXPORT ####
gsea_aki <- gsea_aki_tables %>%
    mutate(db = "Hinze", .before = 1)
names(gsea_aki$gsea_flextables) <- gsea_aki$design

4.4 COLLATE AND SUMMARIZE GSEA RESULTS

GSEA results from each library can now be combined and summarised. Note that pathway annotations need to be interpreted in light of their developement in combination with how GSEA determines their enrichment, and the relevance of enriched pathways must be evaluated in the context of transplantation. For our interpretation, we avail to the pathway annotation as well as core enrichment genes contributing to enrichment via GSEA. We have simplified the interpretation of enriched pathways based on the overlap in biological function of their core enrichement genes.

# JOIN THE GSEA RESULTS ####
gsea <- purrr::reduce(
    list(
        gsea_go,
        gsea_msigdb,
        gsea_wiki,
        gsea_dose,
        gsea_kegg,
        gsea_reactome,
        gsea_aki
    ), bind_rows
) %>%
    mutate(keep = map_dbl(gsea_tables, nrow)) %>%
    dplyr::filter(keep > 0) %>%
    dplyr::select(db, design, genes_gsea, gsea_tables) %>%
    mutate(gsea_tables = map(gsea_tables, mutate, Description = Description %>% as.character()))


# DEFINE INTERPREATION BASED ON PATHWAYS ####
pathway_interpretation <- gsea %>%
    dplyr::select(design, gsea_tables) %>%
    unnest(everything()) %>%
    mutate(
        interpretation = case_when(
            ID %in% c(
                "GO:0006952",
                "GO:0009605",
                "GO:0009607",
                "GO:0043207",
                "GO:0044419",
                "GO:0051707",
                "GO:0006955",
                "GO:0002376",
                "DOID:2914",
                "DOID:77",
                "DOID:1579",
                "R-HSA-168256"
            ) & design == "Baseline_vs_Week24" ~ "immune response",
            ID %in% c(
                "R-HSA-168256",
                "R-HSA-168249",
                "R-HSA-1643685",
                "R-HSA-5663205"
            ) & design == "Baseline_vs_Week52" ~ "immune-related response to injury",
            ID %in% c("AKI", "R-HSA-109582") ~ "response to injury",
            ID %in% c(
                "R-HSA-9006934",
                "R-HSA-597592",
                "R-HSA-5653656",
                "R-HSA-392499"
            ) ~ "response to injury",
            ID %in% c(
                "R-HSA-2262752",
                "R-HSA-8953897"
            ) ~ "response to injury",
        )
    ) %>%
    dplyr::select(design, ID, Description, interpretation)


# WRANGLE THE GSEA TABLES ####
gsea_tables <- gsea %>%
    mutate(
        gsea_tables = pmap(
            list(design, gsea_tables),
            function(design, gsea_tables) {
                gsea_tables %>%
                    mutate(
                        design = design,
                        Description = Description %>% as.character()
                    ) %>%
                    left_join(pathway_interpretation, by = c("ID", "Description", "design")) %>%
                    dplyr::select(ID, Description, interpretation, setSize, NES, p.adjust, core_enrichment) %>%
                    distinct(ID, setSize, interpretation, .keep_all = TRUE)
            }
        )
    ) %>%
    unnest(gsea_tables) %>%
    arrange(p.adjust) %>%
    dplyr::filter(p.adjust < 0.001) %>%
    nest(.by = design, gsea_table = c(-design, -genes_gsea)) %>%
    left_join(gsea %>% dplyr::select(design, genes_gsea) %>% distinct(design, .keep_all = TRUE),
        by = "design"
    )


# MAKE FLEXTABLES ####
gsea_flextables <- gsea_tables %>%
    mutate(
        flextable = pmap(
            list(design, gsea_table),
            function(design, gsea_table) {
                cellWidths <- c(2, 3, 7, 4, 1.5, 1.5, 1.5, 13) # for individual tables up or down
                if (design == "Baseline_vs_Week24") {
                    title <- paste("Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week24 (by FDR)", sep = "")
                } else if (design == "Week24_vs_Week52") {
                    title <- paste("Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between week24 and week52 (by FDR)", sep = "")
                } else if (design == "Baseline_vs_Week52") {
                    title <- paste("Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week52 (by FDR)", sep = "")
                }
                gsea_table %>%
                    mutate(
                        library = db,
                        pathway = Description,
                        "core enrichment genes" = core_enrichment %>% str_replace_all("/", ", "),
                        "n genes" = setSize,
                        NES = NES %>% round(2),
                        FDR = case_when(
                            p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
                        ),
                    ) %>%
                    dplyr::select(library, ID, pathway, interpretation, "n genes", NES, FDR, "core enrichment genes") %>%
                    flextable::flextable() %>%
                    flextable::add_header_row(top = TRUE, values = rep(title, ncol_keys(.))) %>%
                    flextable::merge_v(part = "header") %>%
                    flextable::merge_h(part = "header") %>%
                    flextable::border_remove() %>%
                    flextable::border(border = officer::fp_border(), part = "all") %>%
                    flextable::padding(padding = 0) %>%
                    flextable::padding(j = "core enrichment genes", padding.left = 5) %>%
                    flextable::align(align = "center", part = "all") %>%
                    flextable::align(j = "core enrichment genes", align = "left", part = "body") %>%
                    flextable::font(fontname = "Arial", part = "all") %>%
                    flextable::fontsize(size = 10, part = "body") %>%
                    flextable::fontsize(size = 12, part = "header") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::width(width = cellWidths, unit = "cm")
            }
        )
    )

Summary of all enriched pathways at week 24 (compared to baseline)

# PRINT FLETABLES TO POWERPOINT ####
gsea_flextables %>%
    dplyr::filter(design == "Baseline_vs_Week24") %>%
    pull(flextable) %>%
    pluck(1) 

Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week24 (by FDR)
library	ID	pathway	interpretation	n genes	NES	FDR	core enrichment genes
GO	GO:0006952	defense response	immune response	20	-3.07	3e-05	PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B
GO	GO:0009605	response to external stimulus	immune response	20	-3.07	3e-05	PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B
GO	GO:0009607	response to biotic stimulus	immune response	20	-3.07	3e-05	PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B
GO	GO:0043207	response to external biotic stimulus	immune response	20	-3.07	3e-05	PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B
GO	GO:0044419	biological process involved in interspecies interaction between organisms	immune response	20	-3.07	3e-05	PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B
GO	GO:0051707	response to other organism	immune response	20	-3.07	3e-05	PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B
GO	GO:0006955	immune response	immune response	23	-2.93	5e-05	PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, TRAV12-2, FCGR3A, IDO1, GBP1, GZMB, SH2D1B
GO	GO:0002376	immune system process	immune response	24	-2.75	0.0002	PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, TRAV12-2, FCGR3A, IDO1, GBP1, GZMB, SH2D1B
Reactome	R-HSA-168256	Immune System	immune response	14	-2.51	0.0005	PIK3AP1, PTPRC, TRIM4, KLRF1, LYZ, FCGR1A, KLRD1, GNLY, IL18RAP, CXCL10, FCGR3A, GBP1, SH2D1B
DOSE	DOID:2914	immune system disease	immune response	10	-2.56	0.0007	PRF1, FCGR1A, CXCL9, IL18RAP, CXCL10, CXCL11, FCGR3A, IDO1, GZMB
DOSE	DOID:77	gastrointestinal system disease	immune response	15	-2.59	0.0007	PTPRC, EPB41L3, LYZ, PRF1, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, CXCL11, FCGR3A, IDO1, GZMB
DOSE	DOID:1579	respiratory system disease	immune response	10	-2.56	0.0007	PRF1, FCGR1A, CXCL9, IL18RAP, CXCL10, CXCL11, FCGR3A, IDO1, GZMB

Summary of all enriched pathways at week 52 (compared to baseline)

gsea_flextables %>%
    dplyr::filter(design == "Baseline_vs_Week52") %>%
    pull(flextable) %>%
    pluck(1) 

Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week52 (by FDR)
library	ID	pathway	interpretation	n genes	NES	FDR	core enrichment genes
Reactome	R-HSA-168249	Innate Immune System	immune-related response to injury	20	-3.18	6e-06	LRRFIP1, C1QB, JUN, ATP6V1A, NCF2, CFI, EP300, FCGR1A, CAPZA1, CAPZA2, ADAM10, HSPA1A, STAT6, CALM2, SDCBP, PAFAH1B2, ACTB, RHOA
Reactome	R-HSA-168256	Immune System	immune-related response to injury	35	-3.08	2e-05	LRRFIP1, C1QB, CANX, JUN, ATP6V1A, RNF213, NCF2, PIK3AP1, VIM, CFI, EP300, FCGR1A, CAPZA1, CAPZA2, ADAM10, HSPA1A, IFNGR1, STAT6, CALM2, SDCBP, COL1A2, PAFAH1B2, ACTB, COL1A1, RHOA
Reactome	R-HSA-1643685	Disease	immune-related response to injury	32	-2.92	2e-05	LRRFIP1, PPFIBP1, CANX, JUN, RNF213, CD163, PIK3AP1, TPM3, SIN3A, EP300, FCGR1A, RAN, ADAM10, HSPA1A, IFNGR1, CALM2, CEBPB, SBSPON, COL1A2, ITGB3, ACTB, COL1A1, HNRNPK
Hinze	AKI	Cellular response to acute kidney injury	response to injury	39	-2.75	4e-05	SPTBN1, NRP1, KLF6, HLA-DRB1, DOCK11, DUSP1, CHSY1, LDHA, CDC42SE2, TPM4, PTBP3, ARHGAP29, SYNE2, SERP1, CORO1C, LRRFIP1, PPFIBP1, CANX, JUN, SCOC, RNF213, SPP1, PIK3AP1, VIM, SVIL, TPM3, RAN, ANXA5, ADAM10, HSPA1A, CEBPB, SPARC, ACTB, HNRNPK, SLFN5, RHOA
Reactome	R-HSA-9006934	Signaling by Receptor Tyrosine Kinases	response to injury	19	-2.59	0.0002	ATP6V1A, SPP1, NCF2, SIN3A, EP300, ADAM10, STAT6, CALM2, COL1A2, SPARC, ITGB3, ACTB, COL1A1, RHOA
Reactome	R-HSA-597592	Post-translational protein modification	response to injury	17	-2.68	0.0002	CANX, SPP1, PRMT3, PPP6R3, SIN3A, EP300, CAPZA1, CAPZA2, ADAM10, CALM2, SBSPON, ACTB, HNRNPK, RHOA
Reactome	R-HSA-109582	Hemostasis	response to injury	17	-2.54	0.0005	SIN3A, CAPZA1, CAPZA2, ANXA5, CALM2, COL1A2, SPARC, ITGB3, ACTB, COL1A1, RHOA
Reactome	R-HSA-392499	Metabolism of proteins	response to injury	19	-2.42	0.0005	PRMT3, PPP6R3, SIN3A, EP300, CAPZA1, CAPZA2, ADAM10, CALM2, SBSPON, ACTB, HNRNPK, RHOA
Reactome	R-HSA-5653656	Vesicle-mediated transport	response to injury	15	-2.47	0.0008	SCOC, AP4B1, CD163, PPP6R3, CAPZA1, CAPZA2, CALM2, COL1A2, SPARC, PAFAH1B2, ACTB, COL1A1